#include "merge.hpp"


void merge_usage()
{
    std::cerr << "Usage: comb merge blob-file1 blob-file2 ... blob-fileN" << std::endl;
}

// TODO: fix to handle merge of sequence space files.

void do_merge(int argc, char **argv)
{
    if(argc <= 2 || strcmp(argv[1], "help") == 0)
    {
        merge_usage();
        exit(1);
    }

    std::ifstream f(argv[1], std::ios_base::binary | std::ios_base::in);
    if(!f)
        error("Could not open file: " + std::string(argv[1]));
    unsigned int ftype, fversion;
    boost::iostreams::filtering_istream in;
    in.push(boost::iostreams::gzip_decompressor());
    in.push(f);
    boost::archive::binary_iarchive ia(in);
    ia >> ftype;
    ia >> fversion;

    map_str_vafl_t subs, totl;
    map_str_vafl_t subs2, totl2;
    mapping_type maptype, maptype2;
    
    if(ftype == SNP_CALL_COL_FILE_TYPE)
    {
        color_call_vector_t cv1, cv2;
        bool diploid, dip2;

        if( fversion != SNP_CALL_COL_FILE_VERSION)
            error("bad file version for snp call archive.");
        ia >> diploid;
        ia >> cv1;
        f.close();
        for(int i = 2; i<argc; i++)
        {
            read_col_snp_call_file(argv[i], dip2, cv2);
            if(diploid != dip2)
                error("tried to merge diploid and haploid");
    
            unsigned int k,j;
            for(k = 0; k < cv1.size(); k++)
            {
                // TODO: check that snp k corresponds kn both cases.
                double tot = 0.0;
                for(j = 0; j < cv1[k].probs.size(); j++)
                {
                    cv1[k].probs[j] *= cv2[k].probs[j];
                    tot += cv1[k].probs[j];
                }
        
                for(j = 0; j < cv1[k].probs.size(); j++)
                    cv1[k].probs[j] /= tot;
            }
        }

        write_col_snp_call(std::cout, diploid, cv1);
    }
    else if(ftype == SNP_CALL_SEQ_FILE_TYPE)
    {
        if( fversion != SNP_CALL_SEQ_FILE_VERSION )
            error("bad file version for snp call archive.");

        seq_call_vector_t cv1, cv2;
        bool diploid, dip2;

        ia >> diploid;
        ia >> cv1;
        f.close();
        for(int i = 2; i < argc; i++)
        {
            read_seq_snp_call_file(argv[i], dip2, cv2);
            if(diploid != dip2)
                error("tried to merge diploid and haploid.");

            unsigned int k,j;
            for(k = 0; k < cv1.size(); k++)
            {
                double tot = 0.0;
                for(j = 0; j < cv1[k].probs.size(); j++)
                {
                    cv1[k].probs[j] *= cv2[k].probs[j];
                    tot += cv1[k].probs[j];
                }

                for(j = 0; j < cv1[k].probs.size(); j++)
                    cv1[k].probs[j] /= tot;
            }
        }
        
        write_seq_snp_call(std::cout, diploid, cv1);
    }
    else if(ftype == SUBS_TOTL_FILE_TYPE)
    {
        if( fversion != SUBS_TOTL_FILE_VERSION)
            error("bad file version for subs_totl archive.");

        ia >> maptype;
        ia >> subs;
        ia >> totl;

        f.close();

        for(int i = 2; i < argc; i++)
        {
            read_map_count_file(argv[i], maptype2, subs2, totl2);

            if(maptype != maptype2)
                error("Maptype files are not compatible. bad mapping_type. (mixing color and sequence reads?)");

            std::map<std::string, std::valarray<float> >::iterator mi;
    
            for(mi = subs.begin(); mi != subs.end(); ++mi)
            {
                std::string s = mi->first;

                if(subs.find(s) == subs.end())
                    error("files have different chromosomes!");

                subs[s] += subs2[s];
                totl[s] += totl2[s];
            }
        }

        write_map_count(std::cout, maptype, subs, totl);
    }
    else
    {
        error("Unknown merge file type");
    }
}

